By Lauren Parker, Community Legal Services (4/4/2019)
This report presents data on local and statewide Housing Choice Voucher (HCV) usage to better understand where source of income discrimination may be occurring and support related advocacy efforts.
This report was built using open-source programming, tools and datasets, including R, QGIS, GitHub.
All datasets used in this report are publicly-available and from administrative sources.
# Import Housing Choice Voucher Data
hcv <- read.csv("data/HUDPicture_2017_HCV.csv", stringsAsFactors = FALSE) %>%
filter(Id2 != "42XXX999999") %>%
mutate(tract = Id2) %>%
mutate(county = str_sub(tract, 1, 5)) %>%
select(tract, county, hcv_sub_units)
tenure <- read.csv("data/ACS_17_5YR_B25003_with_ann_TENURE.csv", stringsAsFactors = FALSE, colClasses=c("Id2"="character")) %>%
mutate(tract = Id2) %>%
select(-Id, -Id2)
hcv <- hcv %>%
left_join(tenure, by = "tract") %>%
mutate(hcv_rate = ifelse(tothh < 10, NA, 100 * hcv_sub_units/tothh))
hcv.county <- hcv %>%
group_by(county) %>%
summarise(hcv_sub_units = sum(hcv_sub_units, na.rm = TRUE),
tothh = sum(tothh, na.rm = TRUE)) %>%
as.data.frame() %>%
mutate(hcv_rate = 100 * hcv_sub_units/tothh)
# Import Census Data
hisp <- read.csv("data/ACS_17_5YR_B03002_with_ann_HISP.csv", stringsAsFactors = FALSE)
pov <- read.csv("data/ACS_17_5YR_S1701_with_ann_POV.csv", stringsAsFactors = FALSE)
hisp <- hisp %>%
mutate(tract = as.character(Id2)) %>%
mutate(tract_poc = totpop - nhisp_wh) %>%
mutate(tract_pct_poc = 100*(totpop - nhisp_wh)/totpop) %>%
select(tract, totpop, tract_poc, tract_pct_poc)
pov <- pov %>%
mutate(tract = as.character(Id2)) %>%
mutate(tract_pct_pov = 100 * (belpov/totpovstatus)) %>%
select(tract, totpovstatus, belpov, tract_pct_pov)
# Import Tract Shapefile for Philadelphia
tracts <- st_read("shps/tl_2018_42_tract.shp", layer = "tl_2018_42_tract", stringsAsFactors = FALSE) %>%
mutate(tract = GEOID) %>%
select(tract, geometry)
# Import County Shapefile
counties <- st_read("shps/PA_Counties.shp", layer = "PA_Counties", stringsAsFactors = FALSE) %>%
mutate(state = str_sub(GEOID, 1, 2)) %>%
filter(state == "42") %>%
mutate(county = GEOID) %>%
select(county, geometry, NAME)
# Join data to tracts and counties
tracts <- left_join(tracts, hcv, by = "tract")
tracts[, 3:7][is.na(tracts[, 3:7])] <- 0 # Replace NAs with 0s
tracts <- tracts %>%
mutate(county = str_sub(tract, 1, 5)) %>%
left_join(hisp, by = "tract") %>%
left_join(pov, by = "tract") %>%
mutate(recap = ifelse(tract_pct_poc>=50 & tract_pct_pov >= 40, "RECAP", "NOT RECAP")) # RECAP definition by HUD: https://data.world/hud/recap
counties <- left_join(counties, hcv.county, by = "county")
# Import HUD Small Area FMRs and Zip Codes
safmr <- read.csv("data/hud_phila_small_fmrs_fy2019.csv", stringsAsFactors = FALSE)
fmr <- read.csv("data/hud_phila_fmrs_fy2015.csv", stringsAsFactors = FALSE) #https://www.huduser.gov/portal/datasets/fmr/fmrs/FY2015_code/2015summary.odn
zcta <- st_read("shps/Phila_ZCTA_WGS84.shp", layer = "Phila_ZCTA_WGS84", stringsAsFactors = FALSE)
# Import Median Gross Rent Data
rent <- read.csv("data/ACS_17_5YR_B25064_with_ann_RENT.csv", stringsAsFactors = FALSE) %>%
mutate(high_moe = ifelse(moe_median_rent/median_rent > 0.40, 1, 0)) %>%
mutate(tract = as.character(GEOID))
This section uses 2017 data from HUD to show the spatial distribution of households with vouchers, both state-wide and locally in Philadelphia.
Here is a table of voucher households by county:
county.tbl <- counties
st_geometry(county.tbl) <- NULL
county.tbl <- county.tbl %>%
mutate(hcv = round(hcv_sub_units),
hh = round(tothh),
hcv_rate = round(hcv_rate,2)) %>%
select(NAME, hcv, hh, hcv_rate)
datatable(county.tbl, rownames = FALSE, colnames = c("County", "HCV", "Households", "HCV per 100 Households"),
options = list(order = list(list(1, 'desc'))))
The following is a dot density map of voucher holders across Pennsylvania. Each dot represents 10 households with vouchers and is randomly placed within it’s respective Census tract (the dots do not represent the exact location of households). Click on each county to see the number and percentage of voucher holders.
dots <- suppressMessages(st_sample(tracts, size = round(tracts$hcv_sub_units/10), type = "random")) # each dot = 10 hcv units
popup <- paste0("County: ", counties$NAME, "<br>", "Voucher Rate: ", round(counties$hcv_rate, 2), "<br>", "Vouchers: ", counties$hcv_sub_units)
m <- leaflet(options = leafletOptions(preferCanvas = TRUE)) %>%
addProviderTiles("Esri.WorldGrayCanvas", options = providerTileOptions(
updateWhenZooming = FALSE, # map won't update tiles until zoom is done
updateWhenIdle = TRUE # map won't load new tiles when panning
)) %>%
addCircleMarkers(data = dots, weight = 0, color = "#18BC9B", radius = 3) %>%
addPolygons(data = counties, fill = FALSE, stroke = TRUE, color = "#626468", weight = 1, popup = ~popup) %>%
addLegend("bottomright", colors= "#18BC9B", labels="1 dot = 10 HCVs", title="Vouchers in PA Counties")
m
The following map shows the estimated percent of voucher holders by neighborhood in Philadelphia. Click on individual neighborhoods to see the number and percent of voucher holders in that area. Neighborhood estimates were produced from tract-level estimates by applying area-weights.
# Calculate HCVs by Neighborhood
tract.neigh <- read.dbf("shps/Tracts_Neigh_Int.dbf", as.is = TRUE) %>%
mutate(tract = GEOID) %>%
mutate(neigh = MAPNAME) %>%
mutate(int_fract = AreaInt_ft/Area_ft) %>%
select(tract, neigh, Area_ft, AreaInt_ft, int_fract) %>%
left_join(hcv, by = "tract") %>%
mutate(hcv_weighted = hcv_sub_units * int_fract) %>%
mutate(hh_weighted = tothh * int_fract)
neigh.hcv <- tract.neigh %>%
group_by(neigh) %>%
summarise(
hcv = sum(hcv_weighted, na.rm = TRUE),
hh = sum(hh_weighted, na.rm = TRUE)
) %>%
as.data.frame() %>%
mutate(hcv_rate = ifelse(hh < 10, NA, 100*hcv/hh))
neigh <- st_read("shps/Neighborhoods_Phila_WGS84.shp", layer = "Neighborhoods_Phila_WGS84", stringsAsFactors = FALSE, quiet = TRUE) %>%
mutate(neigh = MAPNAME) %>%
select(neigh, geometry) %>%
left_join(neigh.hcv, by = "neigh")
bins <- quantile(neigh$hcv_rate, seq(0,1, by = 1/5), na.rm = TRUE)
pal <- colorBin("YlOrRd", domain = neigh$hcv_rate, bins = bins)
popup <- paste0("Neighborhood: ", neigh$neigh, "<br>", "Voucher Rate: ", round(neigh$hcv_rate, 2), "<br>", "Vouchers: ", round(neigh$hcv, 0))
m2 <- leaflet(neigh) %>%
setView(lng = -75.140012, lat = 39.999786, zoom = 11) %>%
addProviderTiles("Esri.WorldGrayCanvas") %>% #Stamen.TonerBackground
addPolygons(color = "#ffffff", weight = 1, smoothFactor = 0.5, opacity = 0.8, fillOpacity = 0.5,
fillColor = ~pal(hcv_rate), popup = ~popup,
highlightOptions = highlightOptions(color = "#444444", weight = 2, bringToFront = TRUE)) %>%
addLegend(pal = pal, values = ~hcv_rate, opacity = 0.7, title = "Voucher Holders </br> per 100 Households", position = "bottomright")
m2
neigh.tbl <- neigh
st_geometry(neigh.tbl) <- NULL
neigh.tbl <- neigh.tbl %>%
mutate(hcv = round(hcv),
hh = round(hh),
hcv_rate = round(hcv_rate,2))
datatable(neigh.tbl, rownames = FALSE, colnames = c("Neighborhood", "HCV", "Households", "HCV per 100 Households"),
options = list(order = list(list(3, 'desc'))))
HUD established the definition of Racially and Ethnically Concentrated Areas of Poverty (RECAP) as Census tracts with:
See: https://data.world/hud/recap for more information.
This table presents the number and percent of tenants with vouchers in RECAP tracts vs. non-RECAP tracts.
recap <- tracts %>%
group_by(county, recap) %>%
summarise(
total_households = sum(tothh, na.rm = TRUE),
total_hcv = sum(hcv_sub_units, na.rm = TRUE),
ave_hcv_rate = mean(hcv_rate, na.rm = TRUE),
ave_pct_poc = mean(tract_pct_poc, na.rm = TRUE)
) %>%
as.data.frame() %>%
left_join(counties, by = "county") %>%
filter(is.na(recap)==FALSE) %>%
select(county, NAME, recap, total_hcv) %>%
spread(recap, total_hcv) %>% # convert to wide
as.data.frame() %>%
rename(not_recap = `NOT RECAP`) %>%
rename(recap = RECAP) %>%
mutate(pct_recap = round(100*(recap/(recap+not_recap)),2))
datatable(recap, rownames = FALSE, colnames = c("County ID", "County", "HCVs not in RECAPs", "HCVs in RECAPs", "% HCVs in RECAPs"),
options = list(order = list(list(4, 'desc'))))
*Results for Philadelphia
phl.tracts <- tracts %>%
filter(county == "42101")
phl.recap <- phl.tracts %>%
filter(recap == "RECAP")
bins <- quantile(phl.tracts$hcv_rate, seq(0,1, by = 1/5), na.rm = TRUE)
pal <- colorBin("YlOrRd", domain = phl.tracts$hcv_rate, bins = bins)
popup <- paste0("Voucher Rate: ", round(phl.tracts$hcv_rate, 2), "<br>", "Vouchers: ", round(phl.tracts$hcv_sub_units, 0))
m3 <- leaflet(phl.tracts) %>%
setView(lng = -75.140012, lat = 39.999786, zoom = 11) %>%
addProviderTiles("Esri.WorldGrayCanvas", group = "Base") %>%
addPolygons(group = "Voucher Rate",
color = "#ffffff",
weight = 1, opacity = 0.5,
fillOpacity = 0.5, fillColor = ~pal(hcv_rate),
popup = ~popup,
highlightOptions = highlightOptions(color = "#444444", weight = 2)) %>%
addPolygons(data = phl.recap,
group = "RECAP",
color = "#3b485e", weight = 3,
fill = "#ffffff", fillOpacity = 0.2) %>%
addLegend(pal = pal, values = ~hcv_rate, opacity = 0.7, title = "Voucher Holders </br> per 100 Households", position = "bottomright") %>%
addLegend(data = phl.recap, colors = "#3b485e", labels = "RECAP Tracts", position = "bottomright") %>%
addLayersControl(
baseGroups = "Base",
overlayGroups = c("Voucher Rate", "RECAP"),
options = layersControlOptions(collapsed = FALSE)
)
m3
In late 2016, HUD began using small area Fair Market Rents (FMRs) in certain metropolitan areas to set payment standards for Housing Choice Voucher recipients. In the Philadelphia-Camden-Wilmington metro area, zip codes are used as small area FMRs.
Small area FMRs were an attempt to correct for the error in using county or metro-wide Fair Market Rent estimates. In theory, small area FMRs would be closer to the on-the-ground rent reality; however, zip codes are still large geographies that do not match perfectly to more granular neighborhood-level rental prices. For example, in gentrifying areas, small area FMRs may underestimate rents, further disincentivizing landlords from renting to tenants with vouchers.
The following map shows the difference between 2019 small area FMRs for a 2-bedroom unit and the 2015 metro-wide FMR for a 2-bedroom unit ($1,156). A positive number reflects areas
safmr <- safmr %>%
mutate(efficiency = gsub("\\$", "", efficiency)) %>%
mutate(br_1 = gsub("\\$", "", br_1)) %>%
mutate(br_2 = gsub("\\$", "", br_2)) %>%
mutate(br_3 = gsub("\\$", "", br_3)) %>%
mutate(br_4 = gsub("\\$", "", br_4)) %>%
mutate(efficiency = gsub(",", "", efficiency)) %>%
mutate(br_1 = gsub(",", "", br_1)) %>%
mutate(br_2 = gsub(",", "", br_2)) %>%
mutate(br_3 = gsub(",", "", br_3)) %>%
mutate(br_4 = gsub(",", "", br_4)) %>%
mutate(efficiency = as.integer(efficiency)) %>%
mutate(br_1 = as.integer(br_1)) %>%
mutate(br_2 = as.integer(br_2)) %>%
mutate(br_3 = as.integer(br_3)) %>%
mutate(br_4 = as.integer(br_4)) %>%
mutate(br2_diff = br_2 - 1156)
zcta.fmr <- zcta %>%
mutate(zip = as.integer(GEOID10)) %>%
select(zip, geometry) %>%
left_join(safmr, by = "zip")
bins <- c(-166, -75, -5, 5, 200, 644)
pal <- colorBin("PuOr", domain = zcta.fmr$br2_diff, bins = bins)
popup <- paste0("Difference $: ", zcta.fmr$br2_diff)
m4 <- leaflet(zcta.fmr) %>%
setView(lng = -75.140012, lat = 39.999786, zoom = 11) %>%
addProviderTiles("Esri.WorldGrayCanvas") %>%
addPolygons(color = "#ffffff", weight = 1, smoothFactor = 0.5, opacity = 0.8, fillOpacity = 0.5,
fillColor = ~pal(br2_diff), popup = ~popup,
highlightOptions = highlightOptions(color = "#444444", weight = 2)) %>%
addLegend(pal = pal, values = ~br2_diff, opacity = 0.7, title = "FMR Difference", position = "bottomright")
m4
# lit <- read.csv("SOI Literature Summary.csv", stringsAsFactors = FALSE)
# datatable(lit, rownames = FALSE)
Based on NYC analysis by the Commision on Human Rights and Mayor’s Office of Data Analytics.
NYC Methodology:
Methodology:
phl.tracts.rent <- phl.tracts %>%
left_join(rent, by = "tract") %>%
filter(high_moe != 1)
bins <- quantile(phl.tracts.rent$median_rent, seq(0,1, by = 1/5), na.rm = TRUE)
pal <- colorBin("YlOrRd", domain = phl.tracts.rent$median_rent, bins = bins)
popup <- paste0("Rent: ", phl.tracts.rent$median_rent)
m6 <- leaflet(phl.tracts.rent) %>%
setView(lng = -75.140012, lat = 39.999786, zoom = 11) %>%
addProviderTiles("Esri.WorldGrayCanvas") %>%
addPolygons(color = "#ffffff", weight = 1, smoothFactor = 0.5, opacity = 0.8, fillOpacity = 0.5,
fillColor = ~pal(median_rent), popup = ~popup,
highlightOptions = highlightOptions(color = "#444444", weight = 2)) %>%
addLegend(pal = pal, values = ~median_rent, opacity = 0.7, title = "Median Rent by Tract", position = "bottomright")
## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'
m6